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Abstract 

The study of environmental problems usually requires the description of variables with differ- 
ent nature and the assessment of relations between them. In this work, an algorithm for flexible 
estimation of the joint density for a circular-linear variable is proposed. The method is applied 
for exploring the relation between wind direction and SO2 concentration in a monitoring station 
close to a power plant located in Galicia (NW-Spain), in order to compare the effectiveness of 
precautionary measures for pollutants reduction in two different years. 

Keywords: Circular distributions, Circular kernel estimation, Circular-linear data, Copula. 

1 Introduction 

Air pollution studies require the investigation of relationships between emission sources and pol- 
lutants concentration in nearby sites. In addition, the effectiveness of environmental policies with 
the aim of pollution reduction should be checked, at least in a descriptive way, in order to assess 
whether the implemented precautionary measurements had worked out. 

Different statistical methods have been considered for the study of the relation between wind di- 
rection and pollutants concentration, both for exploratory and for inferential analysis, taking into 
account that wind direction is a circular variable requiring a proper statistical treatment. Wind 
potential assessment using descriptive methods and spectral analysis has been carried out by Shih 
(2008). In addition, wind direction has been proved to play a significant role in detection of emission 
sources (see Chen et al. (2011)) and air quality studies (see Bayraktar et al. (2010)), although the 
wind direction is not treated as a circular variable, but discretized as a factor. Jammalamadaka 
and Lund (2006), consider regression models for the pollutants concentration (linear response) over 
the wind direction (circular explanatory variable), constructing the regression function in terms of 
the sine and cosine components of the circular variable. Recently, Deschepper et al. (2008) intro- 
duce a graphical diagnostic tool, jointly with a test, for fit assessment in parametric circular-linear 
regression, illustrating the technique in an air quality environmental study. 

From a more technical perspective, Johnson and Wehrly (1978) and Wehrly and Johnson (1979) 
present a method for obtaining joint circular-linear and circular-circular densities with specified 
marginals, respectively. Fernandez-Duran (2004) introduces a new family of circular distributions 
based on nonnegative trigonometric sums, and this idea is used in Fernandez-Duran (2007) in the 
construction of circular-linear densities, adapting the proposal of Johnson and Wehrly (1978). The 
introduction of nonnegative trigonometric sums for modelling the circular distributions involved in 
Johnson and Wehrly (1978) formulation allows for more flexible models, that may present skewness 
or multimodality, features that cannot be reflected through the von Mises distribution (the classical 
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model for circular variables). However, the flexibility claimed in this proposal can be also obtained 
by a completely nonparametric approach. 

In this work, a procedure for modelling the relation between a circular and a linear variable is pro- 
posed. The relation is specified by a circular-linear density which, represented in terms of copulas, 
can be estimated nonparametrically. The estimation algorithm can be also adapted to a semipara- 
metric framework, when an underlying model for the marginal distributions can be imposed, and 
it also comprises the classical Johnson and Wehrly model. It also enables the construction of an 
estimation framework without imposing an underlying parametric model. The copula approach 
presents some computational advantages, in order to carry out a simulation study. 




Figure 1: Locations of monitoring station (circle) and power plant in Galicia (NW-Spain). Location 
of station B2: 7° 44' 10" W, 43° 32' 05" N. Power plant location: 7° 51' 45" W, 43° 26' 26" N. 

The practical aim of this work is to explore the relation between wind incidence direction (wind 
blowing from this direction) and SO2 levels in a monitoring station close to a power plant located 
in Galicia (NW-Spain). The monitoring station and the thermal power plant locations are shown 
in Figure 1. In the power plant, energy was usually produced from the combustion of local coal, 
which also generates pollutants as sulphur dioxide (SO2). In order to reduce the emission of SO2 to 
the atmosphere, and to comply with European regulations, coal with less sulphur content has been 
used since 2005. In addition, the power plant changed the energy production system by settling a 
combined process of coal and gas burning in 2008. These measures were aimed to reduce the SO2 
emissions. 

The analized data corresponds to SO2 and wind incidence direction measurements taken during 
January 2004 and January 2011, with one minute frequency. The monitoring station B2 (see Figure 
1) is located in a wind farm 13.4 kilometres in the NE direction with respect to the power plant. In 
Figure 2, rose diagrams for wind direction are shown, including also the average SO2 concentration 
for each wind direction sector. It can be seen that higher values of SO2 are shown in 2004. Note also 
that there are two dominant wind incidence directions, specifically, blowing from SW and from NE. 
In addition, note that the average SO2 values for winds blowing from SW are remarkably high for 
2004, which is explained by the position of the monitoring station with respect to the power plant. 



2 



However, the relation between wind direction and SO2 levels is not clear from these representations, 
and the dependency (or lack of dependency) between them should be investigated. 



January 2004 January 201 1 




Figure 2: Rose diagrams for wind direction in station B2 for 2004 (left panel) and 2011 (right panel), 
with average SO2 concentration. 

This work is organized as follows. Section 2 provides some background on circular and circular- 
linear random variables and a brief review on copula methods. The algorithm for estimating a 
circular-linear density is detailed and discussed in Section 3. The finite sample properties of the 
algorithm, combining parametric and nonparametric approaches, are illustrated by a simulation 
study. The completely nonparametric version of this algorithm is applied for analyzing the relation 
between wind direction and SO2 concentration in Section 4. Some final comments are given in 
Section 5. 

2 Background 

The main goal of this work is to describe the relation between wind direction and SO2 concentra- 
tion in a monitoring station close to a power plant at two different time moments, before and after 
precautionary measurements to reduce SO2 emissions were applied. Bearing in mind the different 
nature of the variables and noticing that measurements from wind direction are angles, some back- 
ground on circular and circular-linear random variables is introduced. This methodology will be 
needed in order to describe the wind direction itself and the joint relation between the two variables. 
For that purpose, the Johnson and Wehrly (1978) family of circular-linear distributions will be in- 
troduced. A general procedure, based on the copula representation of a density, allows for a more 
flexible estimation framework. The goal of this copula representation is twofold: firstly, the classical 
Johnson and Wehrly family can be written in such a way, just involving univariate (circular and 
linear) densities; secondly, with the copula representation, flexible circular-linear relations beyond 
this specific model are also possible. 

2.1 Some circular and circular— linear distributions 

Denote by a circular random variable with support in the unit circle Sr. A circular distribution 
\E r (-) for assigns a probability to each direction (cos(0), sin(0)) in the plane M 2 , characterized 
by the angle 6 G [0, 2tt) (see Mardia and Jupp (2000) for a survey on statistical methods for 
circular data). The von Mises distribution is the analogue of the normal distribution in circular 
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random variables. This family of distributions, usually denoted by vM(fi, k), is characterized by two 
parameters: /J, E [0, 2ir), the circular mean and k > 0, a circular concentration parameter around \x. 
The corresponding density function is given by 

<p vM {P;»,K) = -J— e «^-f), 6 E [0,2vr), (1) 

being io( K ) = 5^ e KCOSUJ duj the modified Bessel function of first kind and order zero. The von 
Mises cumulative distribution function, considering the zero angle as the starting point, is defined as: 



*„m(0;M)«) = / <PvM(u;n,K)dLj, 9e[0,2n). 
Jo 



The uniform circular distribution, 



M«) = ^, ^e[0,2vr), (2) 



is obtained as a particular case of the von Mises family, for k = 0. Circular density estimation 
can be performed by parametric methods, such as Maximum Likelihood, or using nonparametric 
techniques, based on kernel approaches, as will be seen later. 

In order to explain the relation between a circular and a linear random variable (in our example, 
wind direction and SO2 concentration), the construction of a joint circular-linear density will be 
considered. A circular-linear random variable (0, X) is supported on the cylinder S 1 x I or in 
a subset of it and a circular-linear density for (Q,X), namely p(-,-)> must satisfy a periodicity 
condition in the circular argument, that is: 

p{9, x) = p{6 + 2-rrk, x), 6 e [0, 2vr), x G M, k G Z, 

as well as the usual assumptions on taking nonnegative values and integrating one. Johnson and 
Wehrly (1978) propose a method for obtaining circular-linear densities with specified marginals. 
Denote by ip(-) and /(■) the circular and linear marginal densities, respectively, and by *(■) and 
F(-) their corresponding distribution functions. Let also g(-) be another circular density. Then: 

p(6, x) = 2irg [2tt (tf (0) - F(x))] <p(0)f(x), (3) 

is a density for a circular-linear distribution for a random variable (Q,X), with specified marginal 
densities ip(-) and /(•) (see Johnson and Wehrly (1978), Theorem 5). Circular-linear densities with 
specified marginals can be also obtained considering the sum of the marginal distributions in the ar- 
gument of the joining density in (3). Circular-linear densities may include von Mises and Gaussian 
marginals, but the dependence between them will be specified by the joining density g(-). In fact, 
the independence model corresponds to a constant (2tt)^ 1 joining density. From a data sample of 
(G, X), assuming that the joint density can be represented as in (3), an estimator of p(-, •) could be 
obtained by the estimations of the marginals and the joining density. Wehrly and Johnson (1979) 
prove that the construction of circular-circular distributions (that is, distributions on the torus) 
can be done similarly to (3), just considering prespecified circular marginal distributions. Note 
that (3) is just a construction method and not a characterization of circular-linear densities. In 
addition, there are no available testing procedures for checking if a certain dataset follows such a dis- 
tribution. Hence, a more general approach for circular-linear density construction would be helpful. 



In the next section, some background on copulas will be introduced, allowing for a more flexible 
procedure for obtaining circular-linear densities with specified marginals, where the representation 
in (3) fits as a particular case. With this proposal, a fully nonparametric estimation procedure can 
be applied. 
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2.2 Some notes on copulas 

Copula functions are multivariate distributions with uniform marginals (see Nelsen (2006) for a 
complete review on copulas). One of the main results in copula theory is Sklar's theorem, which, 
in the bivariate case, states that if F(-, •) is a joint distribution function with marginals -Fi(-) and 
i^O) then there exists a copula C(-, ■) such that: 

F(x,y) = C(F 1 (x),F 2 (y)), x,y eWL. (4) 

If F±(-) and F 2 (-) are continuous distributions, then C(-, •) is unique. Conversely, if C(-, •) is a copula 
and Fi(-) and F 2 (-) are distribution functions, then F(-,-) defined by (4) is a bivariate distribution 
with marginals F±(-) and F 2 (-). 

If the marginal random variables are absolutely continuous, Sklar's result can be interpreted in 
terms of the corresponding densities. Denoting by c(-, •) the copula density, the bivariate density of 
F(-,-) in (4) can be written as 

f(x,y) = c(F 1 (x),F 2 (y))f 1 (x)f 2 (y), x,yeR. 

As pointed out by Nazemi and Elshorbagy (2012), copula modelling provides a simple but powerful 
way for describing the interdependence between environmental variables. In our particular setting, 
the nature of the variables is different, being the variables of interest linear (SO2 concentration) 
and circular (wind direction). Circular-linear copulas, that will be denoted by Cq^x {.'■>'), can be 
also defined taking into account the characteristics of the circular marginal, satisfying ce ; x(0, v) = 
CQ,xO-,v), Vf G [0, 1], where ce,x( - , ■) is the corresponding circular-linear copula density. Hence, a 
circular-linear density with marginals ip(-) and /(■) is given by: 

p(0, x) = c ,x W), F(x))<p(0)f(x), e G [0, 2vr), x G R. (5) 

Note that Johnson and Wehrly proposal can be seen as a particular case of (5). For a certain joining 
density g(-) in (3), the corresponding copula density is given by: 

c ,x (*(9),F(x)) = 2tt 9 [2vr (tf (0) ± F(x))} , (6) 

where the sign ± refers to the possibility of considering the sum or the difference of the marginal 
distributions in the argument of g(-), as it was previously mentioned. 

A circular-circular density, following Wehrly and Johnson (1979) result, can be done similarly 
just considering circular marginals for (6, $7) and guaranteeing that the copula density satisfies 
ce,n(0,i>) = cs t n(l,v) and CQ t n(u,0) = cq^(u,1), Vu,v G [0,1]. For simplicity, the copula density 
in (6) will be denoted as J&W copula (see Figure 3, left panel). 

The representation of a circular-linear density in (5) enables the construction of new families of 
circular-linear distributions. From Corollary 3.2.5. in Nelsen (2006), a new family of circular-linear 
copulas with quadratic section (QS copula) in the linear component can be constructed. The copula 
densities in this family are given by 

c e,x( u ' v ) = 1 + 2vracos(27ru)(l - 2v), |ct| < (2tt) _1 . (7) 

The QS copula family is parametrized by a, which accounts for the deviation from the independence 
copula corresponding to a = 0. Figure 3 (middle plot) shows the wavy surface corresponding to 
a = (27r)~ 1 . The position of the three modes in the density, centred along u = 0, u = 0.5 and u = 1, 
as well as their concentration, is controlled by the value of a. 
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A possible way to derive new copulas is through mixtures of other copulas (see Nelsen (2006)). 
Thus, for any copula £(•, •), the mixture 

C0,x(u, v) = - (c(u, v) + 5(1 —u,v) + c(u, 1 — v) + 5(1 — u, 1 — f )) , (8) 

leads to a new copula satisfying ce,x(0,u) = ce,x(l> ^)> Vv G [0,1]. This construction also satisfies 
C0,x(w, 0) = ce,x(^> 1)) Vw £ [0, 1], and provides a circular-circular copula (which is also circular- 
linear). A parametrized copula density Cq X (-,-) can be obtained considering, for example, the 
parametric Frank copula: 

- , x a(l - e^V-^"^) 

((1 - e~ a ) - (1 - e- Q «)(l - e"^)) 2 

The mixture copula (8) will be referred to as the reflected copula of c(-, •). The parameter a also 
measures the deviation from independence, which is a limit case as a tends to zero. The copula 
density surface can be seen in Figure 3 (right plot). In this example, the copula density surface 
shows five modes concentrated in the corners and the middle point of the unit square, and the 
peakedness of the modes increases as a grows. 



These three families will be considered in the simulation study. It should be also mentioned that the 
copula representation poses some computational advantages in order to reproduce by simulation data 
samples from circular-linear distributions (see Section 3.1). Finally, although this work is focused 
on the circular-linear case, some comments will be also made about circular-circular distributions. 




Figure 3: Copula density surfaces. Left panel: J&W copula with von Mises joining density with 
parameters \i = 7r and k = 2. Middle panel: QS-copula with a = (27T)" 1 . Right panel: reflected 
Frank copula with a = 10. 



3 Estimation algorithm 

Denote by Xi)}f =1 a random sample of data for the circular-linear random variables (Q,X) 
and consider the copula representation for p(-, •) in (5). In this joint circular-linear density model, 
three density functions must be estimated: the marginal densities (£>(•) and /(•) (and also the cor- 
responding distributions) and the copula density ce,x( - , •)• A new natural procedure for estimating 
p(-, ■) is given in the following algorithm: 

Estimation algorithm 

Step 1. Obtain estimators for the marginal densities </?(•), /(•) and the corresponding marginal 
distributions #(•), F(-). 
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Step 2. Compute an artificial sample { (*(6j) , F(Aj)) } n =1 and estimate the copula density ce,x( - > ■)■ 
Step 3. Obtain the circular-linear density estimator as p(6,x) = ce,x(*(#), F{x))(p{0)f{x). 

The estimation of the marginal densities in Step 1 can be done by parametric methods or by non- 
parametric procedures. For instance, a parametric estimator for /(•) (respectively, for F{-)) can 
be obtained by Maximum Likelihood. In the circular case, that is, for obtaining </?(■), Maximum 
Likelihood approaches are also possible (see Jammalamadaka and SenGupta (2001), Chapter 4). 
These estimators are consistent, although restricted to parametric families such as the von Mises 
distribution or mixtures of von Mises. 

Nonparametric kernel density estimation for linear random variables was introduced by Parzen 
and Rosenblatt (see Wand and Jones (1995) for references on kernel density estimation) and the 
properties of this estimator have been well studied in the statistical literature. Consider {Aj}™ =1 a 
random sample of a linear variable X with density /(■). The kernel density estimator of /(■) in a 
point x G K is given by 

1=1 v ' 

where K(-) is a kernel function (usually a symmetric and unimodal density) and h is the bandwidth 
parameter. One of the crucial problems in kernel density estimation is the bandwidth choice. There 
exist several alternatives for obtaining a global bandwidth minimizing a certain error criterion, usu- 
ally the Mean Integrated Squared Error, such as the rule-of-thumb, least-squares cross-validatory 
procedures (see Wand and Jones (1995)) or other plug-in rules, like the one proposed by Sheather 
and Jones (1991). 

Hall et al. (1987) introduced a nonparametric kernel density estimator for directional data in the 
(/-dimensional sphere E> q . For the circular case (q = 1), denoting by O a random variable with 
density </?(•), the circular kernel density estimation from a sample {6j}" =1 is given by 

, n n 

W) = — #e[0,2vr), (10) 
n i=i 

where L(-) is the circular kernel, v is the circular bandwidth and cq{v) is a constant such that <~p u {') 
is a density. Some differences should be noted in contrast to the linear kernel density estimator 
in (9). First, the kernel function L(-) must be a rapidly varying function, such as the exponential 
(see Hall et al. (1987)). Secondly, the behaviour of v is opposite to h: in linear kernel density 
estimation, small values of the bandwidth h produce undersmoothed estimators (small values of 
v oversmooth the density), whereas large values of h give oversmoothed curves (large values of u 
produce undersmoothing) . See Hall et al. (1987) for a detailed description of the estimator and its 
properties. 

As in the linear case, bandwidth selection is also an issue in circular kernel density estimation. Al- 
though in the linear case it is a well-studied problem, for circular density estimation there are still 
some open questions. Hall et al. (1987) propose selecting the smoothing parameter by maximum 
likelihood cross-validation. There are other recent proposals, such as the automatic bandwidth se- 
lection method introduced by Taylor (2008), but based on his results none of the selectors proposed 
seems to show a superior behaviour. 

Although for the marginal distributions it may be reasonable to assume a parametric model, it is not 
that clear for the copula function, regarding the dependence structure between the marginals. Hence, 
in a general situation, the copula estimation in Step 2 should be carried out by a nonparametric 
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procedure that will be explained below. However, for the Johnson and Wehrly density in (3), the 
copula density ce,x( - , •) is linked with a joining circular density g(-) in (6) and this circular density 
can be estimated in the same way as the marginal circular density. Note that, in this family, all 
the estimators involved in the algorithm are obtained in a strictly univariate way, which simplifies 
their computation. 
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Figure 4: Illustration for data reflection for copula density estimation. The central plot corresponds 
with the original data quantiles. Left panel: mirror reflection from Gijbels and Mielniczuk (1990). 
Right panel: circular-mirror reflection for estimator (11). 



Nonparametric copula density estimation can be also done by kernel methods, as proposed by Gi- 
jbels and Mielniczuk (1990). The proposed estimator is similar to the classical bivariate kernel 
density estimator, with a product kernel and with a mirror image data modification. This mirror 
image (see Figure 4, left panel) consists in reflecting the data with respect to all edges and corners 
of the unit square, in order to reduce the edge effect. In our particular case of circular-linear copula 
densities, the reflection must be done accounting for the circular nature of the first component, as 
shown in Figure 4, right panel. 

The copula density kernel estimator can be defined as: 

n 9 

c & , x (u,v) =T^ B (tt-^ef),.-^)) , (11) 
1=1 1=1 

where and F(-) are the marginal kernel distribution functions, Kb(u, v) = \B\~ l K{B~ 1 / 2 (u, v)) 
is a bivariate rescaled kernel (see Ruppert and Wand (1994) for notation) and B is the bandwidth 
matrix (|_B| denotes the determinant). A plug-in rule for selecting the bandwidth matrix B was 
proposed by Duong and Hazelton (2003). The reflected data sample in each square I = 1,. . . ,9, 
namely { (*(6? } ) , P{xf } )) }™ =1 is obtained by rows in the 3x3 plot (Figure 4, right panel) as 
follows: (u — 1, — v), (u, — v), (u + 1, — v) (bottom row); (u — 1, v), (u, v), (u + 1, v) (middle row); 
(u — 1, 2 — v), (u, 2 — v), (u + 1, 2 — v) (top row). 

As pointed out by Omelka et al. (2009), the estimator proposed by Gijbels and Mielniczuk (1990) in 
the linear case still suffers from corner bias, which could be corrected by bandwidth shrinking. This 
issue is not so important in the circular-linear setting, given the periodicity condition. It should 
be also noticed that, with this construction, the estimator obtained in Step 3, with the proposed 
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reflection, is guaranteed to be a density as long as the distribution and density estimators satisfy 
that F'(-) = /(•) and *'(•) = $(■). 

The estimation algorithm can be extended for circular-circular densities, just with suitable (cir- 
cular) density estimators for the marginals and a slight modification of the data reflection for the 
copula estimation. Specifically, reflection in 3 x 3 scheme as the one shown in Figure 4, the mid- 
dle one corresponding to the original circular-circular data quantiles will be done as follows: (u — 
l,v), (u, v), (u+1, v) (bottom row); {u— 1, v), (u, v), (u+1, v) (middle row); (u— 1, v), (u, v), (it+1, v) 
(top row). 

3.1 Some simulation results 

In order to check the performance of the estimation algorithm for circular-linear densities, the fol- 
lowing scenarios are reproduced. Examples 1 and 2 were proposed by Johnson and Wehrly (1978). 
Example 3 corresponds to the QS-copula family given by (7) and Example 4 to the reflected Frank 
copula constructed in (8). 

Example 1 (J&W copula with circular uniform and Normal marginal distributions). Let (fu(') 
denote the circular uniform density (2) and <f>{-) the standard Normal density ($(•) the standard 
Normal distribution). Take the joining density g{ ) = <Pvm(-', P-, k)- A circular-linear density with 
margins ip(-) and </>(■) is given by 

Pi(9,x) = — exp [kcos(# - 2ir$(x) - //)] <pu(0)(f>(x). 

Example 2 (J&W copula with von Mises and Normal marginal distributions). Consider g(-) = 
<Pvm('', fJ>' , k') an d density marginals </>(■) and <£Vm(s f^2, 1^2) (with corresponding distributions <&(•) 
and ^ v m(-', 1^2), respectively). A joint circular-linear density is given by 

P2(0,x) = - 1 exp [k'cos (2tt(^ vM (0;H2,K2) - ®(x)) ~ A*')] ¥>«m(0; H2, K 2 )<t>(x). 

Example 3 (QS-copula with von Mises and Normal marginal distributions). Take ^VfcfOj K 3) 
and (p(-) as circular and linear marginals. The circular-linear density with copula (7) and a = 
(27r) -1 is given by 

<PvM(0',IJ>3,K3)<f>(x)- 

Example 4 (Reflected Frank copula with von Mises and Normal marginal distributions). Take 
<Pvm{;\ A*4j ^4) and </>(■) as circular and linear marginals. The circular-linear density with reflected 
Frank copula (a = 10) is obtained by the mixture construction (8): 

Pi(6, x) = Cq X (*„m(0; A*4, «4), $(x)) (Pvm(0; /M, k^{x). 

As commented in Section 2, the formulation of the joint circular-linear density in terms of copulas 
simplifies the simulation of random samples. The general idea is to split the joint distribution 
P(-, •) by Sklar's Theorem in a copula Cq : x(-, ■) and marginals and F(-). Simulating a sample 
from (U, V), uniform random variables with copula Ce,x( - 5 ■)> and applying the marginal quantiles 
transformations to get (^~ l (U), F _1 (V r )), we obtain a sample from distribution P(-, •). 



p 3 (9,x) 



1 + — cos(27r*„ M (0; Us, «3))(1 - 2$(x)) 

Z7T 



9 



Figure 5: Density surfaces for the simulation study. Top left: Example 1 with /x = ir and k = 2. 
Top right: Example 2 with // = tt, k! = 5, [i<i = tt/2 and K2 = 2. Bottom left: Example 3 with 
a = (27r) _1 , /i3 = 7r/2 and K3 = 0.5. Bottom right: Example 4 with a = 10, \i± = 7r/2 and K4 = 0.5. 

The simulation of (U,V) values from the copula Ce,x( - > - ) can be performed by the conditional 
method for simulating multivariate distributions (see Johnson (1987)). The conditional distribution 
of V given U = u, denoted by C u (-), can be expressed as 



where the first equality is an immediately property of copulas. So, for simulating random samples 
for the examples, or more generally, for simulating random samples of circular-linear random vari- 
ables with density (5), we may proceed with the following steps: 

Simulation algorithm: 

Step 1 Simulate U ~ 1), W ~ M(0, 1), where U(0, 1) stands for the standard uniform distri- 




(12) 



bution. 



Step 2 Compute V = C' 1 



(W). 
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Step 3 Obtain 6 = * _1 (?7), X = F' 1 ^). 



In Step 1, two independent and uniformly distributed random variables are simulated. The condi- 
tional simulation method for obtaining (U, V) from the circular-linear copula Ce,x(-, ■) is performed 
in Step 2. Finally, quantile transformations from the marginals are applied, obtaining a sample from 
(@,X) following the joint density (3). For Examples 1 and 2, the conditional distribution C u (-) in 
(12) is related to a von Mises vM(ll, k) and in Step 2, for each u, V = (2-7r)~ 1 (W; 2iru — ll, k). 

For Example 3, for each u, V = a+1 ^aW ^ w ^ a _ 2iracos (2iru). Example 4 is 

simulated from the mixture of copulas (8), using that for the Frank copula, Step 2 becomes 



The estimation algorithm proposed will be applied to estimate the densities in the examples. Dif- 
ferent versions of the estimation algorithm can be implemented, considering parametric and non- 
parametric estimation methods. In addition, for the Johnson and Wehrly densities (Examples 1 and 
2), the estimation of the circular-linear density can be approached by representation (3), in terms 
of a circular joining density, or by the more general representation (5). Summarizing, the following 
variants of the estimation algorithm will be presented: 

I Johnson and Wehrly, parametric [JWP]: based on representation (3), marginals as well as 
joining density are parametrically estimated. The results will be used as a benchmark for the 
Johnson and Wehrly models (Examples 1 and 2). 

II Johnson and Wehrly, semiparametric [JWSP]: based on representation (3), marginals are 
estimated parametrically and a nonparametric kernel method is used for the joining density. 

III Johnson and Wehrly, nonparametric [JWNP]: based on representation (3), marginals and 
joining density are estimated by kernel methods. 

IV Copula, semiparametric [CSP]: based on representation (5), parametric estimation is consid- 
ered for marginals. The copula density is estimated by kernel methods. 

V Copula, nonparametric [CNP]: based on representation (5), marginals and copula density are 
estimated by kernel methods. 

In the parametric case, density estimators have been obtained by Maximum Likelihood, specifying 
the von Mises family for the circular distributions and the Normal family for the linear marginal. 
Nonparametric estimation has been carried out using kernel methods. The kernel density estimator 
in (9), with Gaussian kernel and Sheather- Jones bandwidth, has been used for obtaining /(•). For 
(£>(•) and g(-), the circular kernel density (10) has been implemented, with exponential kernel and 
likelihood cross-validatory bandwidth. In the semiparametric approaches (parametric marginals 
and nonparametric joining density or copula), Maximum Likelihood has been used for obtaining 
(£>(•) and /(•)• The circular kernel estimator (10) has been considered for g(-). The copula density 
kernel estimator (11) has been used for CQ t x(-,-), with bivariate Gaussian kernel and restricted 
bandwidth matrix. Specifically, following the procedure of Duong and Hazelton (2003), full band- 
width matrices have been tried, but with non significant differences in the values of the principal 
diagonal. Hence, a restricted bandwidth with two smoothing parameter values is used, considering 
the same element in the principal diagonal and a second element in the secondary diagonal, regard- 
ing for the kernel orientation. 

In order to check the performance of the procedure for estimating circular-linear densities, the Mean 
Integrated Square Error (MISE) is considered: 
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Johnson and Wehrly Copula-based Relative efficiency 



Ti 
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JVV or JVVlu v^O-T ^lN.r 
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1UUU 


0.0053 0.0074 0.0139 
0.0027 0.0042 0.0085 
0.0005 0.0011 0.0024 

U.UUUO U.UUUU U.UUl^i 


0.0142 0.0185 
0.0092 0.0118 
0.0033 0.0039 
n nn99 n nn9zt 


0.7209 0.3831 0.3773 0.2888 
0.6324 0.3123 0.2879 0.2252 
0.5042 0.2319 0.1647 0.1400 
u.^oiu u.iyui u.izou u.iuyo 


Ex. 2 50 
100 
500 
1000 


0.0406 0.0467 0.0831 
0.0209 0.0252 0.0560 
0.0043 0.0061 0.0180 
0.0021 0.0034 0.0106 


0.0750 0.0824 
0.0485 0.0538 
0.0157 0.0175 
0.0094 0.0105 


0.8691 0.4884 0.5415 0.4927 
0.8304 0.3730 0.4313 0.3887 
0.7029 0.2387 0.2739 0.2457 
0.6271 0.1987 0.2241 0.2016 



Table 1: MISE for estimating the circular-linear density in Examples 1 and 2. Estimation methods: 
Johnson and Wehrly (JWP, JWSP, JWNP) and copula-based (CSP, CNP). Relative efficiencies for 
JWSP, JWNP, CSP and CNP estimations with respect to JWP. 



The MISE is approximated by Monte Carlo simulations, taking 1000 replicates. Four sample sizes 
have been used: n = 50, n = 100, n = 500 and n = 1000. In the first example, the set-up parame- 
ters are fx = ir and k = 2. For the second example, we take // = tt, k' = 5, (12 = tt/2 and K2 = 2. 
Both in Examples 3 and 4, the parameters in the von Mises marginal were set to /i^ = ^4 = tt/2 
and K3 = K4 = 0.5, and the linear marginal is a standard normal. 

In Figure 5, surface plots for the example densities are shown, the top row corresponding to the 
Johnson and Wehrly family (Example 1 and Example 2). Simulation results for these examples can 
be seen in Table 1. For all the alternatives of the algorithm, the MISE is reduced when increasing 
the sample size. Example 2 presents higher values for the MISE, and it is due to the estimation of 
a more complex structure in the circular marginal density (circular uniform in Example 1 and von 
Mises in Example 2). In both models, the estimation methods providing the information about the 
Johnson and Wehrly structure (JWP, JWSP, JWNP) (that is, based on representation (3)) work 
better, as expected. Nevertheless, the copula based approaches, CSP and CNP, are competitive 
with the JWNP. 

The parametric method JWP presents the lowest MISE values for all sample sizes in both examples, 
so it will be taken as a benchmark for computing the relative efficiencies of the nonparametric and 
semiparametric approaches, both based on the density (3) or on the copula density (5). Relative 
efficiencies are obtained as the ratio between the MISE of the parametric method and the MISE 
of the nonparametric and semiparametric procedures. The relative efficiencies (see Table 1) are 
higher for the semiparametric approach, with better results for Example 2. Boxplots for the ISE for 
sample size n = 500 can be seen in Figure 6, top row. The larger variability of the nonparametric 
methods (JWNP and CNP) can be appreciated for both examples. For Example 2 (see Figure 6, 
top right panel), note that JWNP shows the highest values for the ISE. 

Obviously, the first three proposals (JWP, JWSP, JWNP) make sense for distributions belonging 
to Johnson and Wehrly family. However, since the marginals considered in Examples 3 and 4 also 
belong to the von Mises (circular) and the Normal (linear) classes, as in Examples 1 and 2, JWSP, 
JWNP, CSP and CNP variants of the algorithm will be applied in the last two cases. Results are 
reported in Table 2. MISE values decrease with sample size, showing CSP and CNP a similar 
behaviour. Note also that assuming a representation like (3), increases dramatically the MISE: for 
instance, in Example 3, for n = 500 or n = 1000, the MISE for JWSP or JWNP is four times the 
one provided by CNP. The effect of this misspecification in the copula structure can be clearly seen 
in Figure 6, bottom row. 
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Table 2: MISE for estimating the circular-linear density in Examples 3 and 4. Estimation methods: 
Johnson and Wehrly (JWSP, JWNP) and copula-based (CSP, CNP). 



Example 1 Example 2 




JWP JWSP JWNP CSP CNP JWP JWSP JWNP CSP 



Example 3 Example 4 




JWSP JWNP CSP CNP JWSP JWNP CSP 



Figure 6: ISE for Example 1 (top left), Example 2 (top right), Example 3 (bottom left) and Example 
4 (bottom right) for n = 500 and different estimation procedures. 
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4 Application to wind direction and SO2 concentration 



The goal of this work is to explore the relation between wind incidence direction and SO2 con- 
centration in monitoring station B2 near a power plant (see Figure 1 for location of station B2). 
SO2 is measured in ng/m 3 and wind direction as a counterclockwise angle in [0, 27r). With this 
codification, 0, |, tt and ^ represent east, north, west and south direction, respectively. 

The dataset contains observations recorded minutely in January 2004 and January 2011, but due to 
technical limitations in the measuring device, SO2 is only registered when it is higher than 3[ig/m 3 . 
Concentration values below this threshold are considered as non significant and are recorded as the 
lower detection limit (3/ig/m 3 ). Data have been hourly averaged, resulting 736 observations for 
2004 and 743 for 2011. In order to avoid repeated data, perturbation procedures have been applied 
to both marginals, and will be detailed below. Afterwards, a Box-Cox transformation for the SO2 
concentration with A = —0.84 for 2004 and A = —7.34 for 2011 have been applied. For the sake of 
simplicity, we will refer to these transformed data as SO2 concentration, but note that figures are 
shown in the transformed scale. 

Measurement devices, both for the wind direction and for SO2 concentration, did not present a 
sufficient precision to avoid repeated data, and this problem was inherited also for the hourly 
averages. The appearance of repeated measurements posed a problem in the application of the 
procedure, specifically, in the bandwidth computation. Perturbation in the linear variable, the SO2 
concentration, was carried out following Azzalini (1981). A pseudo-sample of SO2 levels is obtained 
as follows: 

Xi = Xi + bei, 

where Xi denote the observed values, b = 1.3crn _1//3 and e^, i = l,...,n are iid random variables 
from the Epanechnikov kernel in (— \/5, \/5). a is a robust estimator of the variance, which has 
been computed using the standardized interquartile range. Azzalini (1981) shows that this choice 
of b for the data perturbation allows for consistent estimation of the distribution function, getting a 
mean squared error with the same magnitude as the one from the empirical cumulative distribution 
function. 

The problem of repeated measures also occurs for wind direction. In this case, a perturbation 
procedure similar to the linear variable case can be used, just considering the circular variable (the 
angle) in the real line. Then, a pseudo-sample of wind direction was obtained as 

6i = 6i + dei, 

with 0i denoting the wind direction measurements and £j, i = 1, . . . , n were independently generated 
from a von Mises distribution with \i = and k = 1, with d = n -1 / 3 . We have checked by simulations 
that the applied perturbation did not affected the distribution of the data. 

4.1 Exploring S0 2 and wind direction in 2004 and 2011 

The estimation procedure is applied to data from 2004 and 2011 in station B2, considering the 
fully nonparametric approach. Specifically, a nonparametric kernel density estimator for the SO2 
concentration was used, with the plug-in rule bandwidth obtained by Sheather- Jones method 
(h = 1.75 x 10~ 6 for 2011, h = 0.02 for 2004). For the wind direction, circular kernel density 
estimation has been also performed, with likelihood cross-validatory bandwidth (y = 46.41 for 
2011, v = 78.56 for 2004). In Figure 7, the estimation of the circular-linear density surfaces, with 
the corresponding contour plots, is shown. 
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For 2011, two modes in the SW and NE directions (see Figure 7, right column) can be identified, 
both with a similar behaviour and with quite low values of SO2 concentrations. Recall that the 
scale is Box-Cox transformed (for data in the original scale, see Figure 2). A different situation 
occurs in 2004 (see Figure 7, left column) where, in addition to the two modes that appear in 
2011 for low values of SO2, a third mode arises. This mode is related to winds blowing from SW 
(from the thermal power plant) and to SO2 values significantly higher than for 2011. This relation 
suggests that the additional mode represents SO2 pollutants coming from the power plant, and its 
disappearance for 2011 illustrates the effectiveness of the control measures applied during 2005-2008 
to reduce the SO2 emissions from the power plant. 





Figure 7: Circular-linear density estimator for wind direction and SO2 concentrations in monitoring 
station B2. Left column: 2004. Right column: 2011. 



5 Final comments 

A flexible algorithm for estimating circular-linear densities is proposed based on a copula density 
representation. The method provides a completely nonparametric estimator, but it can be modified 
to accommodate the classical Johnson and Wehrly family of circular-linear distributions. In the 
purely nonparametric version of the algorithm, circular and linear kernel density estimators have 
been used for the marginals, although other nonparametric density estimators could be considered. 
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In addition, the extension of the algorithm for circular-circular density estimation is straightforward. 

In our air quality data application, the precision of the measurement devices posed some extra prob- 
lems in the data analysis. The lack of precision resulted in the appearance of repeated values for 
the wind direction, and a data perturbation procedure was needed in order to apply the algorithm. 
The perturbation method proposed has been checked empirically, and it is inspired by the results 
for kernel distribution estimation, but our guess is that similar results could be obtained with just 
perturbing the data by summing errors from a highly concentrated distribution (e.g. a von Mises 
distribution with large k). Nevertheless, data perturbation in the circular setting needs further 
investigation. Another possible problem that may be encountered in practice, for linear variables, 
is censoring, that may be due to detection limits or other phenomena. Under censoring, the obser- 
vation values are only partially known, and suitable estimation procedures for density estimation 
with censored data should be applied. 

Finally, a natural criticism of our analysis of the air quality data is that the measurements were 
not independent of one another. Temporal dependence could be accounted for in the estimation 
procedure by using a proper bandwidth selector, such as a A;-fold cross- validatory bandwidth for 
the linear kernel density estimator. However, there are not such alternatives for circular data (up 
to the author's knowledge) and the study of bandwidth selection rules for circular dependent data 
is beyond the scope of this paper. 

The simulation study and real data analysis has been carried out in R 2.14 (R Development Core 
Team (2011)), using self-programmed code and package circular (Agostinelli and Lund (2011)). 
For the real data analysis, the computing time for 2004 is 30.77 seconds, taking the computation of 
the copula estimator 8.58 seconds. The same procedures take 31.57 seconds for the 2011 data. All 
computations were done on a computer with 1.6 Ghz core. This shows that the computational cost 
of the method is not high and its application is feasible in practice. 
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